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I*  Introduction 

Probabilistic  processing  is  a  method  of  processing  the  output 

data  of  an  array  of  seismometers  with  the  aim  of  detecting  earthquake 

or  underground  nuclear  blast  signals  in  the  presence  of  ambient  seismic 

noise.  This  method  is  based  on  the  assumption  that  the  array  output 

is  Gaussian  with  mean  zero  and  known  covariance  matrix  0^  or  0 g  depending 

on  the  absence  or  presence  of  signal. 

The  decision  regarding  presence  or  absence  of  the  signal  i%  therefore, 

made  equivalent  to  testing  the  hypothesis  that  the  observed  date  is  from 

a  Gaussian  population  with  covariance  matrix  0^  with  the  simple  alternative 

hypothesis  that  the  covariance  matrix  is 

The  univariate  situation  is:  a  single  value  x  is  observed  from  a 

Gaussian  population  with  mean  zero.  It  must  be  decided  whether  x  is 

2  2 

more  representative  of  a  Gaussian  population  with  variance  or  o2  • 

The  two  density  functions  are 


fx(x) 

1 

up  (-x2/ao2) 

/5n  o. 

f200 

1 

«xp  (-x2/aoj) 

ysr  a2 

(1) 
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where  0^  >  0^  .  Clearly  it  is  more  likely  that  x  is  from  f^(x)  if 
I x |  <  c  and  more  likely  that  x  is  from  fg(x)  if  |x|  >  c.  So  the  method 
of  processing  in  this  case  reduces  to  just  squaring  the  observation  x 
and  comparing  the  result  with  a  fixed  constant. 

If  a  vector  x  of  dimension  k  is  observed,  then  the  two  density 
functions  become: 

fx(x)  -  (2n)“k/2  inj"172  exp  (-xT  O"1  x/2) 


f2(x)  «  (2n)’k/2  lOgl-1^2  exp  (-xT  x/2). 


(2) 


These  equations  can  be  compared  to  find  the  set  of  vectors  x  for  which 

f. (x)  >  f 0(x).  Thus  I f 

j.  d 

xT(0*1-  Og1)  x  >  c  -  In  (I  f^l/if^l )  (3) 

then  fg(x)  >  f^x). 

The  test  statistic  x^(f)^  -  O^Jx  must,  therefore,  be  computed  from 
the  data.  The  detailed  calculation  of  the  test  statistic  in  terms  of 
array  data  is  given  in  the  following  section. 

A  more  rigorous  development  of  the  test  statistic  based  on  Bayes 
Theorem  has  been  presented  in  a  preceeding  special  reoort^ 


II.  Calculation  of  the  Test  Statistic 
A.  Notation 

The  data  from  an  array  of  seismometers  can  be  represented  by 
the  matrix  X  where  the  rows  correspond  to  the  seismometers  and  the  columns 


“  '  ““  tU  •'’,“<‘''u*,‘niMl4,"H0*t‘*l*J'*  *<•  (t(H !i u»M>l I li iiliil 
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to  the  time  sampled  date 


*11  x12  *  *  X1N 


*21  x22  * 


''HI  XM2  *MN 

« 

Thus,  x.j  is  the  observation  from  seismometer  i  at  the  time  j. 

The  elements  of  the  matrix  X  can  be  "strung  out"  into  a  single  vector 
y  in  the  following  manner: 

T  ,  T  T  Tv 

y  -  (x1#  Xg,  •  •  *  xN) 

where  x^  ■  (x, ,,  xol,  •  •  •  ,  xu>).  This  vector  y  is  now  the  vector  of 

J  <i  nJ 

observations  from  v*hich  the  test  statistic  (equation  3)  will  be  computed. 

Since  the  mean  of  y  Is  assumed  to  be  zero,  the  covariance  matrix  of  y  is 
just  _ 

oy  *  yyT 

and  the  dimension  of  0  is  HN  by  MN  since  there  are  HN  elements  in  the 

y  _  _ 

vector  y.  Let  r|j.j|+1  *  x.xj  if  j  >  i  and  r|j.j|+  ±  ■  x{xj  if  j  <  i  and  write 


2N 


yy 


x^  X1X3  *  *  *  X1XN 

x2xi  x2x2  x2x3  •  •  •  x2xN 


xiT 

Vi 


X  X 

N  N 


00 


so  that 


0. 


rl  r2  *  *  ’  rN 


2  rl  *  *  ’  N-l 


where  each  submatrix  r(  |S  an  H  by  H  matrix 


A.  8. 
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B. 


|J 

Inverse  of  R 

The  inverse  of  R  cen  be  obtained  from  the  solutions  of  the  two 


systems  of  equations 


where  the  dot  indicates  a  matrix  of  zeros  and 
a  .  (a"'1)*1 
b1  -  0  -  BTAB. 

From  equation  {*?)  the  last  H  columns  of  (RN)“a  are 


,  1  ^n****0  w  liJ 


Now  partition 


Jl— tT— 

c!  r"’1 
L  •  J 


Oenota  by  tha  ijth  M  by  H  aubnmtrix  of  tha  invaraa  of  HN.  From  (6)  and  (7) 


niA 
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the  following  equations  can  be  written: 


ij  Ij  rN-i-fl  w-j+1 

I#  j  ■  1,  .  »  •  f  N-l 

*N  ■  I*  fp  r*1  T  (q) 

*1*1,  j+1  *ij  4|+l'V  j+1  *9' 

j  “  •  •  •»  N“1 

Now  by  subtracting  (8)  from  (9)*  the  basic  iterative  formula  of  equation  (10) 


is  obtained 


jn  ‘  ’ij  *  Cjlx  (10) 


and  the  entire  inverse  matrix  can  be  c*nerated  starting  with  the  first 
M  columns  of  the  inverse  which  is  the  first  column  of  the  s^j's  available 
from  equation  (5). 


C.  The  Quadratic  Processor 

Suppose  the  matrix  X  ■  (x^,  *  .  *  ,  x^)  is  observed.  Then  the 


test  statistic 


P(y,)  -  y,CyJ 


can  be  evaluated  for  i  -  1,  .  .  .  ,  T-N+l  where 

yj  -  (xj,  ...  .  ,  xJH+fml) 

CO*1-  O'1 

A  program  to  perform  this  calculation  has  been  written.  In  addition. 


an  integrated,  souared  multichannel  Wiener  filter  output  is  available  as 
an  option  In  the  quadratic  processor  program  to  give  a  comparison  between  the 
two  detection  methods* 


A.  B. 


D.  Multichannel  Time  Series  Oata  Generation 


It  is  desired  to  generate  multichannel  data  xj  =  (x^.,  x^,...,  x^. ) 


.... 


with  the  covariance  matrix  R  .  Let  ij  i  ■  2,  ...»  N  and  P^  satisfy  the 
system  of  equations  (5).  Suppose  y|  is  a  supply  of  vectors  satisfying 


Ey.  *  0 


Ey,y.  «  i. 


Then  x.  will  be  generated  by  the  following  equation 

x.  -  -<r|J)T  -{rJJ)T  x|  N  +  C,  (n) 

for  i  >  N  where  x^  *  x^  «=  c^,  .  .  .  ,  x^_1  ®  The  m*thod  °f 

determining  the  €{  is 

*,  =  Hy, 

where  HH^  ■  PN«  The  reason  for  this  is  that  €j  must  satisfy 

E€i€i  =  PN* 

Determination  of  an  H  satisfying  the  above  is  accompl ished  by  requiring 
that  H  be  lower  triangular. 

The  x  data  generated  by  Equation  II  has  the  correlation  matrix  R 
of  Equation  5«  This  is  established  by  multiplying  Equation  11  on  the 
right  by  x^_{  and  then  taking  the  expected  value. 


This  equation  was  given  incorrectly  in  Section  IV  of  the  Array 
Research  Semiannual  Report  No.  2. 


i 
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E.  Eigenva’ues 

The  test  statistic  is 

P(x)  «  xT  (o'1  -  O^1  )x  (12) 

where  P^  and  P  are  Toepiitz  matrices.  Let  the  observation  vector  x  be 
transformed  to  z  by  the  NM  by  NM  nonsingular  transformation  x  =  Sz  such 
that 

P(x)  -  zTST(P’1  -  P'3 )  Sz  =  zTDz  (13) 

where  0  is  diagonal.  Then  the  computation  of  the  quadratic  form  is  reduced 
to  filtering  operations,  followed  by  a  sumption  of  the  squared  outputs.  A 
especially  important  transformation,  diagonalizing  both  0^,  0^  simultaneously, 
exists  since  A  *  0^  and  B  «  are  positive  definite.^  This  transformation 
Is  found  by  solving  the  following  generalized  eigenvalue  problem. 

The  solution  for  X's  satisfying 

I A  -  XB|  -  0  (14) 

is  called  the  generalized  eigenvalue  problem.  The  corresponding  set  of 
vectors  x  such  that 

Ax  *  XBx  (15) 

are  the  generalized  eigenvectors.  The  matrix  M  whose  columns  consist  of 
these  eigenvectors  normalized  ro  that  xTBx  -  1  can  then  be  shown  to  satisfy 

MTAM  «  A 

MTBM  -  I  (16) 

where  A  is  a  diagonal  matrix  with  elements  X. 

Thus  if  the  S  of  equation  (13)  is  taken  to  be  H,  then  the  diagonal 
matrix  0  satisfies 

0  »  A  -  I.  (17) 


**  “  >  ’»  3?  1%. 

if*  t  =r  'V 
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.II.  Theoretical  Evaluation 
A.  Quadratic  Processor 

It  Is  deslrabie  to  determine  operating  parameters  for  the.  quadratic 
detection  function 

p(x)  -  XT  (n*1  -  n*1}*.  (18) 


Two  useful  parameters  are  the  false  alarm  rate  a  and  the  failure  to  detect 
rate  8.  These  parameters  are  defined  In  terms  of  the  distribution  of  P(x) 
when  x  Is  from  say  fj,  and  when  x  Is  from  fig,  say  fjj.  The  covariance 
matrix  0^  will  represent  noise  and  0g  will  represent  signal  plus  noise.  The 
parameters  a,  8  can  now  be  expressed  by 


a  »  Prob  /  P(x)  >  c  |  x  c  I  \ 

/  /  x  .  (  <19> 

8  -  Prob  ^  P(x)  <  c  I  x  6  II) 

and  are  the  probabilities  of  the  two  types  of  error  which  can  be  made.  The 
other  situations 

P(x)  <  c  I  x  €  I 

,  (20) 

P(x)  >  c  j  x  €  II 

represent  correct  decisions  by  the  method  of  probabll Is'lc  processing. 

Suppose  the  data  x  Is  transformed  by  the  matrix  H  of  section  E. 


Then 


x  »  Hy 
P(x)  - 


yV  (a;1  - 

yT  (A  -  I)y 


(2!) 
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where  A  Is  a  diagonal  matrix  with  elements  Xj  satisfying  the  generalized 
eigenvalue  equation 


-  0. 


(22) 


It  follows  from  the  definitions  of  0^  and  that  X>1.  Anrther  equivalent 
expression  for  22  Is 

I  0,  c^1  -  XI  |  -  0  (£3) 

so  that 

la,"*1  I  ■  l  °2 1  -  "  X|  •  W 

Thus  the  constant  c  of  Equation  3  can  be  written 


C  ■  In  II  X  -  E  In  X  • 
i  '  I  1 


(25) 


If  we  Insert  ■  N  and  ♦  N  into  Equation  23  It  follows 

that 

i  Slf 1  -  (X  -  V  I  |  -  0.  (26) 

Since  S  and  N  are  positive  definite,  the  roots  of  SN*1  are  positive  and 
therefore  X j>1. 

The  distribution  of  y,  If  x  c  II,  is  multivariate  normal  with  mean 

-IT  -1 

zero  and  covariance  matrix  M  M  •  I»  It  follows  from  Equation  21  that 

P(x)  has  the  same  distribution  as 


* 


fe7) 


. . .  .  ,  - . .  J 
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wnere  ■  Xj  -  1  and  are  Independent  chi-square  random  variables 

each  with  1  degree  of  freedom. 

If  x  €  I,  then  y  is  multivariate  normal  with  covariance  M*^  P.,  M  ^ 
A  Now  transform  y  to  a  by 

y  -  A “1/2  a  (28) 

so  that  covariance  of  a  is  I  and 


so 


P(x)  -  yT  (A  -  l)y  -  «T  A"1/2  (A  -  I)A‘1/2  k 
-  »T  (I  -  A'1), 

that  P(x)  again  has  the  distribution  as  the  expression  (27)  with 


(29) 


^  (30) 

This  result  could  also  be  established  by  defining  the  columns  of 
M*  to  be  the  eigen  vectors  corresponding  to  the  eigenvalues  A'  satisfying 

I  n"1  -  a  •  n’1  |  -  °  , 

where 

M'T  fi^1  M'  -  I 

M,T  n“X  M*  «  A' 
and  If  x  «  H*y  ,  then 

P(x)  -  yT  (I  *  A')y  •  Now,  when  x  c  I,  y  Is  multivariate 
normal  with  mean  zero  and  covariance  matrix  I  so  that  P(x)  Is  distributed 
as  Equation  27  with  -  1  -  A '  but  A •  -  1/A  yielding  the  seme  result 
as  before. 


d 


5^3- 
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The  exact  distribution  of  Equation  27  is  not  known  in  closed  form.  Several 

3 

approximations  have  been  considered  in  the  literature* 

We  wi 1 1  be  content  to  assume  that  Equation  27  is  approximately  distributed 
as  ^  where  6,  v  are  determined  by  setting  the  mean  and  variance  of  6ic2  equal 
the  mean  and  variance  of  Equation  27.  This  results  in 

4  -  a2/ 2p 

(31) 

V  »  2p2/02 

where  P 

p  =  L  t| 

1-1  1 

p  (32) 

CT2  =  2  Z 

1-1  ' 

It  Is  possible  to  arrive  at  the  above  approximation  by  applying  a 

jl 

result  of  Patnalk  ,  which  relates  a  non-central  chi-square  variable  to  a 
constant  multiple  of  a  chi-square 


(33) 


to  each  term  of  Equation  27;  thery  use  the  addition  theorem  for  non-centrai  chi- 


square  values#  and  finally  apply  Patnalk's  result  to  the  summed  non-central 


chi-squares*  The  advantage  of  the  above  approach  Is  in  transmitting  the 
known  qualities  of  the  approximation  (33)  to  the  approximation  given  by 


Equations  31  snd  32* 

The  above  results  can  be  applied  to  the  pair  of  Equations  (19) 


to  obtain 


to  taniwM!  !i  mi  i« ! >«* 
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a  »  Prob  ^ y -2  >  c/6  I  tl|  ■  1  ■  X’-^ 
0  »  Prob  <  c/6  |  T|j  -  X  -  1  y 


(30 


where  v  and  6  are  defined  In  terms  of  the  values  by  Equations  31  end 


10 

• 


Note  that  both  the  approximate  and  true  means  and  variances  are 


given  by 


■ 

Mr  »  ^  (1  *  *j) 

1  1-1  ' 

P 

MtI  -  ^  (*,  “  1) 

^  1-1  1 


r 

Oy  “  2  L  (l  -  X' 

1-1  ' 

P 

0^  -2  I  (X  -  l)2 
J“L  1-1  1 


(35) 


and  that  since  each  term  of  Pj  is  less  than  the  corresponding  term  of 

2  2 

It  follows  that  >  Pj  t  Ojj  >  Oj  as  expected* 

It  Is  not  necessaiy  to  calculate  the  Individual  X's  to  apply  the 
approximation  or  to  obtain  the  parameters  given  by  Equations  35*  From 
section  E, 

•S1 

so  that 

M*1 


M  •  n’1  M  A 

(36) 

n^1  h  -  a 

(37) 
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and 


Also 


tr  A  -  tr  fig  • 

tr  AA  -  tr  M*1  ng  f^1  ng  (l^1  H 


- tr  ^  %  <? 


and  It  follows  that 


1 

I 


S' 


M 


2 


-  tr  °i  "S1  °i  "S1 


(33) 


(39) 


(^0) 


Thus  the  approximation  parameters  v,  6,  and  the  population  means  and 
variances  can  be  written  In  terms  of  the  original  correlation  matrices 
and  their  Inverses. 

B.  Squared  Muiticha.mel  Wiener  Filter 

The  optimum  Wiener  coefficients  to  estimate  the  signal  at  a 
selected  seismometer  are  given  by 


02  W  -  Y  (41) 

where 

*  covariance  matrix  of  signal  plus  noise  (NM  by  NM) 

W  *  optimum  filter  coefficients  (NM  by  1) 

Y  ■  selected  signal  correlations  (NM  by  1). 

The  original  array  data  Is  then  filtered 


N 

-  £ 

1-1 


V' 


‘t  -  I 


WTx 


(42) 
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where 


and  the  resulting  time  trace  is  the  optimum  (mean  square  error) 

estimate  of  the  signal  trace*  It  is  desirable  to  consider  the  output 

2 

squared  ,  as  a  competitor  of  P(x)  for  signal  detection*  We  therefore 

2 

want  th<  istribution  of  y(  when  x^  €  I  and  when  x^ <  II*  Since  y(  is 
a  linear  combination  of  normal  random  variables*  y  Is  normal  with  mean 
zero  and  variance 

o2  -  wT  nx  W.  (43) 

Thus  y2  has  the  same  distribution  as  a2**.2  .  Since  W  -  fl"1  y,  o2  can 
be  written 

O2  -  yt  o'1  nx  n"1  Y  .  (44) 

If  x  c  II*  then  fl  ■  0^  and 

o2  -  yT  Clj1  Y  -  yT  W  (45) 

whereas 

o  -  vT  n*1  ^  nj1  y  (46) 

when  Ml,  The  quant  I ty  y'  W  It  simply  releted  to  th*  mean  square 
error  which  is  output  In  the  multichannel  filter  program  which  determines 
the  filter  coefficients  W* 


-l  6- 

IV.  Discussion  of  Processing  Results 

Correlation  matrices  for  an  annular  ring  noise  model  (3-4  km/sec) 
and  a  solid  disk  (6  km/sec)  signal  model  were  generated  by  an  existing 
program^  and  used  as  input  data  to  the  sequence  of  programs  discussed  in 
Section  II.  The  frequency  spectrum  for  both  signal  and  noi~e  was  left 
white.  The  center  and  three  outer  seismometer  locations  of  the  triangular 
WMSO  array  were  used  thereby  generating  4  by  4  correlation  submatrices 
for  25  lags. 


1 

Geometry  of  the  WMSO  Array  Locations  Used. 

The  10  lag  matrices  for  both  noise  and  signal  plus  noise  were  inverted  and  the 
P,T  matrices  for  each  were  written  on  CPT  to  be  used  in  the  data  generation 
program.  The  random  error  values  are  different  (independently  generated) 
for  each  of  the  6  traces  shown  in  Figures  1  and  2.  The  difference  of  the  2 
inverses  was  also  written  on  CPT  to  be  used  in  the  program  which  evaluates  the 
quadratic  processor.  A  plot  of  the  quadratic  processor  output  as  a  function  of 
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time  is  given  in  figure  1.  The  left  and  right  portions  are  quadratic 
processor  output  for  dat<,  generated  with  the  noise  correlation  matrix 
and  signal  plus  noise  correlation  matrix  respectively.  The  horizontal 
line  separating  the  two  curves  is  drawn  at  the  amplitude  c  of  Section  II. 

The  mean  and  variance  (average  sum  of  squares  of  deviations  about  the 
mean)  calculated  from  the  two  plots  of  figure  1  and  corresponding  theoretical 
values  are 


ml 

') 

51 

m2 

2 

52 


16.0 

\ij  =  17.0 

11.8 

Oj  «=  16.0 

44.4 

Mjj  =  44.1 

217.4 

Ojj  =  227.4 

2 

The  values  of  p  and  o  very  accurately  predict  the  calculated 
2 

voiues  m  and  s  .  These  results  show  that  the  effectiveness  of  the 
quadratic  processor  may  be  conveniently  determined  from  t.e  correlation 
matrices  on  a  theoretical  basis. 


There  are  4  points  of  the  480  noise  points  greater  than  c^0  and  8 
points  of  the  440  signal  plus  noise  points  less  than  c1Q  (ignoring  the 
first  20  points  of  the  second  trace  due  to  an  unfortunate  end  effect). 

The  measured  values  of  a,  8  are  approximately  .006  and  .018  respectively. 
The  theoretical  values  of  a,  B  are  given  by 

a  =  Prob  >  50^  =  .014 

8  =  Prob  <  10}  <  .14  . 

2 

Exact  values  are  unavailable  for  a,  8  due  to  inadequate  %  tables. 
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The  same  models,  using  25  lags,  were  also  evaluated  and  the  results 
are  summarized  in  Table  1.  The  quadratic  processor  output  is  also  plotted 
in  Figure  1. 

The  theoretical  values  for  a,  P  for  25  lags  are  approximately 

a  =  Prob  *q0  >  120  =  .0013 

P  =  Prob  <  20.4  *  .005 

* 

where  the  normal  approximation  has  been  used. 

From  the  evaluation  of  Equations  45  and  46,  it  was  found  that 
for  25  lags  the  squared  four  channel  Wiener  filter  is  distributed  as 

.65  ^  if  x  c  II 

.18  7^  if  X  €  I. 

The  filters  were  determined  and  applied  to  data  generated  from  populations 
I  and  II.  These  results  are  plotted  in  Figure  2.  The  sample  means  of 
I  and  II  are  .72  and  .19* 

No  natural  critical  level  c  is  available  as  for  the  quadratic 
detector.  For  comparison  suppose  we  consider  a  **  B.  Then  the  critical 
level  is  approximately  .15  and  a  =  P  =  .36. 

It  should  be  noted  in  the  above  comparison  that  no  summation 
over  the  squared  output  has  been  made.  The  reason  for  this  is  that 
(25)  •  (4)  of  the  original  data  points  wH  »  produce  one  output  point  by 
either  of  the  two  detection  methods.  If  one  sums  over  a  gate  length  other 
than  1  for  the  Wiener  trace,  then  more  data  is  used  and  it  would  only 
be  fair  to  allow  smoothing  for  the  quadratic  processor. 
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One  argument  against  the  above  is  that  summing  for  one  of  the 
two  methods  may  be  mere  beneficial  than  for  the  other.  The  "whiteness" 
of  the  Wiener  outputs  suggests  smoothing  is  of  greater  importance  to  the 
Wiener  process  than  to  the  quadratic  process.  The  results  of  Figures 
.1  and  2  for  25  lags  were  smoothed  over  a  gate  length  of  10  and  25 
points  and  the  results  shown  in  Figures  3  and  4. 


10  LAGS 


S  +  N 


Figure  1.  Evaluation  of  the  Quadratic  Processor  for  S  =  8  km/ sec  Solid 
Disk  Signal,  N  =  3-4  km/sec  Annular  Ring  Noise,  4  Channels 
(Center  and  Outer  Locations  of  the  WMSO  Array) 


C  =  67.4 


SIGNAL  +  NOISE 


h  SMOOTHED  OVER  10  POINTS 


C*67.4 


SIGNAL  +  NOISE 

■\  SMOOTHED  OVER  25  POINTS 
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Figure  3.  Results  of  Figure  1  Smoothed  Over  10,  25  Points, 
Quadr  .ic  Processor,  25  Lags 
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Figure  4.  Results  of  Figure  2  Smoothed  Over  10,  25  Points, 
Wiener  Filter  Sum  and  Square  Process 
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Table  1.  QUANTITATIVE  RESULTS  FOR  25  LAGS 


I 

II 

m 

43. 56 

336.28 

,2 

97.73 

1570.08 

42.76 

148.90 

o2 

47.20 

1020.17 

101 

35.44  xio‘13 

16.84  X1017 

trace 

trace 

trace  ftj*  Og  O^1 
trace  fig1 


67-45 

248.90 

57.24 

1417.96 


38.08 
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